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Abstract 

Here we present a systematic study of supercoil formation in DNA minicircles under varying 
linking number by using molecular dynamics simulations of a two-bead coarse-grained model. Our 
model is designed with the purpose of simulating long chains without sacrificing the characteristic 
structural properties of the DNA molecule, such as its helicity, backbone directionality and the pres- 
ence of major and minor grooves. The model parameters are extracted directly from full- atomistic 
simulations of DNA oligomers via Boltzmann inversion, therefore our results can be interpreted as 
an extrapolation of those simulations to presently inaccessible chain lengths and simulation times. 
Using this model, we measure the twist /writhe partitioning in DNA minicircles, in particular its 
dependence on the chain length and excess linking number. We observe an asymmetric supercoil- 
ing transition consistent with experiments. Our results suggest that the fraction of the linking 
number absorbed as twist and writhe is nontrivially dependent on chain length and excess linking 
number. Beyond the supercoiling transition, chains of the order of one persistence length carry 
equal amounts of twist and writhe. For longer chains, an increasing fraction of the linking number 
is absorbed by the writhe. 
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A. Introduction 



Conformational features and mechanical properties of DNA in vivo (such as supercoil 
formation, bend/twist rigidity) play an important role in its packing, gene expression, pro- 



tein synthesis. 



protein transport. 



3| etc. Advances in single-molecule probing and 
monitoring techniques in the last decade have provided new opportunities for detailed anal- 
ysis of such mechanical and structural properties. The mechanical response of single DNA 
molecules, {4, 5] the lifetime of denaturation bubbles, 6] and the details of supercoil forma- 



tion 



7|, [s] can now be investigated. These new experimental findings also triggered renewed 



theoretical interest in DNA physics 
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The degree of supercoiling in DNA depends on various tunable parameters, such as the 
linking number, external torque and salt concentration. These parameters essentially mod- 
ify the molecule's relative preference for twisting vs writhing. DNA minicircles provide a 
convenient setting where mechanical properties, such as the asymmetry of the molecule un- 
der positive and negative supercoiling, the possibility of torsion-induced denaturation and 
kink formation can be studied both experimentally |7| and theoretically. 15|, [ly, [l7| Despite 
steady progress in the field in recent decades, a full comprehension of these phenomena is 
still ahead of us. Analytical models have been successful in explaining many qualitative 
aspects of DNA mechanics, however, they are typically too simplistic to capture fine details, 
especially in the short chain limit. Furthermore, their treatment gets difficult in the presence 
of thermal fluctuations and nonlinearities. 

Computer simulations provide a wealth of information on short DNA chains. State of 
the art full-atomistic simulations are capable of microsecond scale simulations of short 
DNA oligomers. [181 Simulations of DNA oligomers have been successfully used to con- 
struct a database 19| analogous to the crystal structure database. Beyond oligomeric DNA 
molecules, full-atomistic simulations have also been performed for DNA minicircles. isl 
However, they are short of providing a through exploration of the conformational space 
even for minicircles as small as 100-300 basepairs (bps). 

In recent years a number of coarse-grained models have been proposed to study the 
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conformational and mechanical features of DNA at intermediate length scales. 
These models are typically generated in an ad-hoc manner with tuning parameters to match 
some of the key features of DNA, such as the pitch, persistence length, melting temperature. 
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and sequence specificity. Drukker&Schatz [20| studied DNA denaturation of sliort B-DNA 
oligomers (10-20 bps) using a two-bead model and no major and minor grooves. Knotts et 



al. [2l| introduced a three-bead model to study the melting dynamics in chains of length ~ 60 
bps. Their model also yields a persistence length of ~ 20 nm compared to the experimentally 
measured value of ~ 50 nm. A more recent model 22| where they study DNA renaturation 



events is also tuned to reproduce the persistence length. Trovato et al. recently proposed 
a single-bead model to study thermal melting that also exhibits major and minor grooves. 
In simulations of 92 and 891 bps long DNA chains, they demonstrated supercoiling and 
denaturation using one sample with positive and negative torsional stress for each length. 
Here, we present our results on the twist/writhe partitioning using a novel coarse-grained 
DNA minicircle model. Using this model we perform a systematic study of the supercoiling 
behavior of DNA, and in particular we investigate the equilibrium amounts of twist and 
writhe accommodated by the chain as a function of the applied torsional stress and the 
chain length. 

Unlike earlier Gb-like approaches, the model parameters are extracted from full-atomistic 
DNA simulations via Boltzmann inversion, with no fitting for structural or mechanical prop- 
erties. With only two beads the model captures most structural details of DNA, such as: 

• the helicity and the pitch, 

• backbone directionality, 

• major and minor groove structure which results in the anisotropic bending rigidity, 

• persistence length. 

The accuracy of our coarse-grained model is mostly limited by the accuracy of the force-field 
used in the full-atomistic DNA simulations. The coarse-graining method employed here 
can be extended in a straightforward manner to include further details, such as basepair 
specificity, hybridization, and explicit charges. However, setting the efficiency of the model 
as the priority for mesoscale simulations, we postpone these extensions to a future study. 
In the next sections, we outline the model and then present our results on the twist/write 
partitioning in DNA minicircles. 
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FIG. 1: A schematic representation of the mapping from the full-atomistic DNA structure (left) 
to a coarse-grained DNA model (right). Note that some of the bonds are omitted in the rightmost 
figure for visual clarity. 

B. Model 

1. Coarse-graining procedure 

The model is composed of two types of "superatoms" P and B per nucleotide, repre- 
senting the collective motion of the backbone phosphate group + sugar (P) and the nucleic 
acid base (B) as depicted in Fig. [H Since a realistic description of the helical structure is 
sought within a minimalist setting the B superatoms are considered generic, with no base 
specificity. The extension to four different bases is straightforward with an approximately 
4-fold increase in the number of model parameters. Above simplification is compensated 
by choosing the superatom positions optimally, with the criterion that the equilibrium dis- 
tributions associated with the degrees of freedom of i?-superatoms have maximal overlap 
when they are calculated for purines {A and G) and pyrimidines (T and C) separately. This 
condition is best met when the Cartesian coordinates of P superatoms are chosen as the 
center of mass of the atoms { 03', P, OIP, 02P, 05', C4', 04', CI', C3', C2' }, while B 
superatoms are placed at the center of mass of the atoms {N9, C8, N7, C5, C6, N3, C4} for 
purines and { Nl, C6, C5, C4, N3, C2} for pyrimidines. 
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FIG. 2: Coarse grained model of DNA molecule based on the super atoms P and B, where the 
former represents the phosphate backbone and the sugar group, and the latter represents the 
nucleic-acid base. The superatoms Pi, Bi from the first strand, and the superatoms P2n-i, B2n-i 
from the second strand form the nucleic-acid basepair which are connected by hydrogen bonds in 
the original system. The intra-strand bonds PiBi, BiPi^i, PjPj+i, BiBi^i are shown by solid lines. 
The inter-strand bonds BiB2n~i and PiP2n~i are shown by dashed lines. (pp^pB ^hpBP 
the dihedral angles defined by Pj, Bi, Pj+i, Pi+i and Bi, Pj+i, Pj+i, Pi+2, respectively. Similarly, 
Oppp and Op^p represent the bond angles defined by Bi^i, Pi, Bi and Pi, Bi, Pi+i. The dihedral 
angle stiffness is explicitly included in the coarse-grained potential, whereas the bond angle stiffness 
is mostly due to the intra-strand PP and BB bonds. 

Fig. [2] shows the coarse-grained model composed of the superatoms P and B and the 
adopted indexing convention for the two strands. 

2. Interactions 

The effective interactions incorporated into the model are shown in Fig. [2l We use four 
bonded and two dihedral potentials that maintain the local single-strand geometry: 
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• harmonic bonds PiBi, BiPi^i, Pj-Pj+i and Bi^iBi that fix the intra-strand superatom 
distances as well as the angles Op^p and Oppp, 

• dihedral potentials associated with the angles 4>pbpb (f^BPBP- 
All four harmonic bond potentials have the form 



where the stiffnesses, Kb, and the equilibrium bond lengths, tq, differ as listed in Table [B 
In particular, the difference between Pj-Bj and BiPi^i bond parameters reflects 5'-3' direc- 
tionality of the molecule. The choice of harmonic bond potentials PiPi+i and Bi^iBi over 
true angular potentials increases computational efficiency without significantly distorting 
the equilibrium distributions, as verified in Fig. [3] for the BPB bond angle. 

The torsional stiffness of the dihedral angles <Ppbpb ^bpbp defined, respectively, by 
the superatoms PiBiPi+iBi^i and Bi^iPiBiPi^i is modeled by the potential. 



where, again, the two stiffness coefficients, Kd, for BPBP and PBPB dihedral angles and 
their equilibrium values, (po, are separately determined from the full- atomistic simulation 
data as described below. 

In addition to the intra-strand interactions above, we define two inter-strand potentials 
that stabilize the double-stranded structure. First is a tabulated potential connecting the 
pairs BiB2n~i and refiects the hydrogen bonding between the nucleic-acid bases A-T and 
G-C (again, base specificity is omitted at this stage). The use of a tabulated potential may 
facilitate monitoring base-pair breaking events in future studies. The second inter-strand 
interaction takes into account the steric hinderence of the base atoms and the electrostatic 
repulsion of the phosphate groups through a repulsive potential between the superatoms 
Pi and P2n~i- This interaction also helps maintain the directional nature of the hydrogen 
bonding between the base pairs, which in our model is displayed by the alignment of the 
superatoms Pi, Bi, B2n-i, and P2n-i- Functional forms of both interactions are discussed in 
the next section. 

Finally, a Lennard- Jones excluded- volume potential between all superatom pairs (except 
PiPi+2) that do not interact via previously defined potentials maintains the self-avoidance 
of the DNA chain. 



(1) 



Vd{4>) = Kd[l - cos{4> - 4>o)] , 



(2) 
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FIG. 3: Potential of mean force curves for intra- and inter-strand interactions obtained by Boltz- 
mann inversion of the relevant distance/angle probability distributions. The target distributions 
obtained from the full-atomistic trajectories (filled black diamonds), the coarse-grained model 
potentials derived from these (solid black lines), and the results of the coarse-grained model sim- 
ulations (open red circles) for the intra- (top) and inter-strand (bottom) degrees of freedom. PB 
bond, PEP angle, and PBPB dihedral angle distributions (not shown) yield similar results. 

3. Determination of the force constants 

The force constants and the equilibrium values for bond and dihedral potentials are 
obtained from the thermal fluctuations of the associated superatoms via Boltzmann inver- 



sion. 



23| The fluctuation data is obtained from molecular dynamics (MD) trajectories of a 



full-atomistic study by Dixit et al. |l9j The full-atom MD data includes room-temperature 
simulations of all possible tetramers of the nucleic-acid bases located at the center of a 15 bps 
long B-DNA oligomers. Since our model does not include base-pair specificity, all tetramer 
data was given equal weight throughout our analysis. 

Boltzmann inversion of the probability distributions for each type of bond length, bond 
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angle, and dihedral angle (obtained from the thermal fluctuations of superatom centers) 
yields potentials of mean force (PMF) which are shown with solid diamond symbols in 
Fig. El 

These PMF curves are used to obtain the force constants for the intra-strand bonded 
interactions by means of harmonic flts. The flts obtained for r^p, and (pBPBP are shown 
with solid-black lines in Fig. [3] and the force constants obtained by these flts are listed in 
Table [H The anharmonic features in PMF curves can be captured using more sophisticated 
potentials, which could be implemented if a more speciflc model is desired. 



Interaction Type 


Equilibrium Position 


Force Constant 


PiBi bond 


ro=5.45 A 


^5:^=7.04 kBT/A^ 


BiPi+i bond 


ro=6.09 A 


Kb=lQ.U kBT/A^ 


PiPi+i bond 


ro=6.14 A 


Kb=20.36 ksT/A^ 


BiBi^i bond 


ro=4.07 A 




PBPB dihedral 


00=3.62 rad 


Kd=25A0 ksT/rad? 


BPBP dihedral 


(/)o=3.51 rad 


Kd=27M kBT/ra(f 



TABLE I: Force constants obtained by fitting the analytical forms in Eqs.([T]) and ([2]) to the Boltz- 
mann inverted distributions shown in Fig[3l 

Also shown in Fig. [3] are the PMF curves for the inter-strand bond distances between 
BiB2n-i and PiP2n~i, where a strong asymmetry is evident in both. The interaction among 
BiB2n~i superatoms stems from the hydrogen bonds as discussed above, and displays an 
equilibrium separation. The BiB2n-i interactions is incorporated into the model via a tab- 
ulated potential, which is obtained by a smooth curve flt to the PMF data. 

The PiP2n-i PMF curve also displays an equilibrium separation. Unlike the BiB2n-i in- 
teraction, here we only model the repulsive part of this interaction via a tabulated potential. 
The attractive part is already captured by the previously discussed potentials, as we will 
discuss below. 
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The excluded volume of the superatoms is represented via repulsive Lennard- Jones inter- 
actions, 

j4[(?^)''-(^)' + 0.25] r<r.,, 
ULAr) = ^ (3) 

measured in units of /c^T. For all B-B and B-P pairs = 5.35 A and r^ut = 6A, whereas 
for all P-P pairs these values are doubled. Superatom pairs that are bonded and all PjPi+2 
pairs are excluded from these Lennard- Jones interactions. The constants and exclusions for 
these Lennard- Jones interactions are chosen such that intra- and inter-strand interactions 
previously defined are not influenced when the Lennard- Jones interactions are switched on. 

The coarse-grained model potential does not include any explicit electrostatic interac- 
tions, therefore strictly speaking this model is suitable for high salt concentrations. 

4- MD simulation 

ESPResSo package 2J] was used for all coarse-grained MD simulations with the Langevin 
thermostat at room temperature. All super-atoms were assumed to have the same mass 
roughly equivalent to 170 atomic mass units. The length and energy units were set to 1 
A and ksT (T=300 K), respectively, from which the unit time (r) can be determined by 
dimensional analysis to be approximately 7 fs. The equations of motion were integrated by 
using the velocity Verlet algorithm with a time-step of 0.015 r. 

We simulated both linear and circular DNA chains. The initial configurations were chosen 
as their respective ground-states obtained through over damped MD simulations. The data 
collection was performed after thermal equilibration, where the required equilibration time 
depended on the measured quantity. For example, the equilibrium distributions for the 
degrees of freedom of the coarse-grained model (Fig. [3]) took approximately 10 CPU minutes 
on an Intel Quad-Core machine using a single processor, whereas thermal averages of the 
writhe and the twist (Fig. [6]) for the longest circular DNA we considered required upto 
64 CPU weeks using 8 processors in parallel. MD trajectories are visualized by VMD 
package. 
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Side view 




FIG. 4: Ground state structure of the model DNA. and top (right) view of the equihbrium structure. 
Key structural features of the DNA molecule, such as the directionality of the backbone and major 
and minor groves are present in the coarse grained model. Note that some of the bonds are omitted 
in the figure for visual clarity. 

5. Equilibrium Properties 

The ground-state structure of the molecule is shown in Fig. HI Despite the fact that 
only two beads are used for each unit, the model successfully captures most of the essential 
features of the DNA molecule. The directionality of the strands is reflected in the differences 
between the force constants for BP vs. PB bonds and PBPB vs. BPBP dihedral angles. 
The double helical structure with major and minor grooves is also captured. The helical 
pitch (the method of calculation is discussed further below) is 11.4 bps in ground state and 
drops to 11.1 bps at room temperature, suggesting an anharmonic twist rigidity which is a 
complex function of the model potentials given above. The reason for the somewhat higher 
helical pitch we find here in comparison with the actual DNA is discussed at the end of this 
section. 

We tested the derived coarse-grained potentials by performing an MD simulation of a 36 
bps long linear chain. The resulting Boltzmann inverted distributions of the bond lengths, 
bond angles and dihedral angles are shown with open red circles in Fig. [31 The distri- 
butions associated with intra-strand interactions nicely match the potentials used, which 
themselves are the best harmonic fits to the corresponding full-atomistic data. The fact 
that the given potentials are recovered from the thermal fluctuations of the superatoms 
suggests that the degrees of freedom used in the coarse-grained Hamiltonian are minimally 
coupled. Also shown in Fig. [3] is the inverted form of the BPB bond angle distribution (top 
row, in the middle). The agreement between the full-atomistic (filled black diamonds) and 
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the coarse-grained (open red symbols) simulation results is observed also for other angular 
potentials (not shown) and justifies the aforementioned use of harmonic bond potentials as 
a replacement for true angular potentials. 

On the other hand, the inter-strand degrees of freedom, namely BiB2n~i and PiP2n-i, 
display some degree of coupling to other bonded interactions. BiB2n~i interaction repre- 
sents the hydrogen bonding within a base-pair. The tabulated potential we used faithfully 
reproduces the equilibrium fiuctuations upto a few ksT (see bottom-left graph in Fig. [3]). 
Higher-energy excitations of the BiB2n~i bond are suppressed beyond the level imposed by 
the derived potential, suggesting a coupling with the other degrees of freedom in this regime. 

PiP2n-i interaction we have used is a purely repulsive interaction, as described above. 
The associated tabulated potential had to be iteratively softened until a good match was 
obtained with the full-atomistic data on the left of the Boltzmann inverted distributions 
(for r < 17.5 A in the bottom-right graph in Fig. [3]). The attractive right-hand-side of 
the effective PiP2n-i potential derived from the superatom fluctuations is solely due to the 
remaining interactions. 

We have looked at the persistence length of our model DNA by means of extended 
simulations of linear molecules with freely fluctuating ends and lengths 114, 228, and 456 
bps. Snapshots from these simulations were used to calculate the persistence length and the 
pitch as explained below. 

Let Qi be the vector joining Pi and Pj+i and C(r) = (ql-qi^r)- Fig. [5]shows C(r) for a 114 
bps long DNA, calculated by a running average along the chain which is further averaged 
over 200 snapshots. The helicity of the molecule results in a sinosoidal function with an 
exponentially decaying amplitude using which both the persistence length and the pitch can 
be measured. The solid line in Fig. [5] represents the flt to C{r) via the following empirical 
function, 

/(r) = e""/'^ [a + (1 - a) cos(27rr/A)] , (4) 

where a is a geometrical constant related to the aspect ratio of the helix and given by 
a = (1/2) — 1/[1 + (A/27ri?)^], A is the helical pitch and R is the radius measured between 
the helix center to the P superatoms. The decay rate of the correlations measured by fltting 
Eq. (jll) in the interval < r < 30 gives the persistence length (Ip) of the model DNA as 96 
bps consistently for chain lengths of 114, 228, 456 bps (data not shown for 228 and 456 bps 
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r{bps) 

FIG. 5: The correlation of the PP vectors as a function of the vector length r and functional fit 
to obtain the persistence length and pitch of the model molecule. 

long chains). At larger chemical lengths (r > 30), finite size effects (short chains) and long 
equilibration times (long chains) limit the accuracy of the data. 

The persistence length of our model DNA is shorter that the experimentally measured 
length of 50 nm's (~ 150 bps). The stiffness of a DNA chain is a result of both the 
bonding/angular interactions that form the local helical structure and the electrostatic self- 
repulsion imposed by the high line charge density. The electrostatic interactions are included 
here only implicitly, through the coarse-grained force-field parameters that best fit the full- 
atomistic simulations performed with explicit electrostatic interactions. This approximation 
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is one possible reason for the shorter persistence length we find, while the imperfection 
(below) of the full-atomistic force-field used by Dixit et al. 19| is another. 

Equation (jlj) simultaneously provides an estimate for the helical pitch at room temper- 
ature as A = 11.1 bps. This value is also slightly higher than that measured 10.5 bps) 
for the B-DNA. An inspection of the full-atomistic data reveals that the same mismatch 
exists for the oligomers simulated by the AMBER parm94 force-field 26| used by Dixit and 



coworkers. 



19| Therefore, the higher pitch is, in fact, a consequence of the imperfect AM- 



BER parm94 force-field rather than a shortcoming of the inverted Boltzmann method used 
here. 

C. Twist vs Writhe in DNA minicircles 



Our main goal in this study is to analyse the twist/writhe partitioning in DNA minicircles. 
In addition to serving as a demonstration of the model's capabilities, this problem is also 
relevant to the denaturation behavior of DNA chains under conserved linking number. 13|, 



Il4 | In this section, we will explain how we measure the response of a circular DNA chain to 
applied torsional stress and present our results for chains of different lengths under varying 



stress levels. 



1. Calculating twist and writhe on a discrete chain 



Twist and writhe reflect two geometrically distinct modes of response of a DNA chain 
to applied torsional stress. Let ri(s) and r2{s) be the two closed curves interpolating B- 
superatoms of each strand, parametrized by the continuous variable s and obtained here by 
the cubic spline method. The centerline of the DNA is given by r(s) = [ri{s) + r2(s)]/2. 
The unit tangent and normal vectors at any point s are 



t{s) = dr{s)/ds 

u{s) = (ri(s) - r2(s))/|ri(s) - r2(s)| 



(5) 



Then, the twist (Tw) of the chain which is a measure of the sum of the successive basepair 



stacking angles is formally given by: 



1 



23] 



ds tis) 



u[s] X 



du{s) 
ds 



(6) 
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The writhe (Wr) is a nonlocal property associated with the torsional stress stored in the 
conformation of the centerline (as in the coiling of the old telephone chords) and can be 
obtained using 

Wr = ^ ([ds([ ds' tis) X tis') ■ ) _ 

For a DNA chain constrained to have a fixed "linking number" (L/c), the number of times one 



chain loops around the other, total twist and total writhe are connected by the relation 28|, 



A variety of met 
discrete chains. 30, 



Lk = Tw + Wr . (8) 
lods have been proposed for calculating the amount of Tw and Wr on 



3l[ | We found that, constructing the curves ri^2{s) by using cubic splines 
and numerically evaluating Eq. ([6]) and Eq. ([7]) is the most accurate approach that guarantees 
validity of Eq. ([8]) at all times. Each snapshot that was used for calculating the average 
values of Tw and Wr plotted in Fig. [6] was checked to satisfy Eq. ([8]) with a percentage error 
< 10~^. Note also that, with the present definition of the centerline as the midpoint of Bi 
and B2n-ii a linear DNA chain (closed at infinity) has a finite writhe density measured as 
Wro/Lk ^ 0.06. This is due to the corkscrew motion of this centerline; a consequence of the 
fact that the inter-strand BiB2n-i bonds do not cross the center of the tube which tightly 
encloses the equilibrium structure. 



2. Twist and writhe under torsional stress 



On circular DNA chains (such as plasmids) the linking number is a topological invariant, 
unaltered by thermal fluctuations {in vivo, topoisomers are employed for this reason). The 
free energy of a circular DNA of length L is minimized when Lk = Lko ^ L/A, where, in a 
strict sense, the equality is attained in the limit L oo and T ^ 0. A weak dependence of 
A on L is possible, but will be ignored here and A is set to the value we found for our linear 
model DNA. 

In this section, we analyze the behavior of DNA minicircles under varying linking number. 
By Eq. ([H]), an excess or a deficiency in the linking number {ALk = Lk — Lk^ ^ 0) modifies 
the equilibrium values of writhe and twist measured in the relaxed state. The partitioning of 
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FIG. 6: The ensemble averaged writhe (black) and twist (red) densities of circular DNA molecules 
for different excess linking numbers cr = [Lk/Lko) — 1. Also shown in green is the sum of the two, 
confirming the agreement with Eq. ([8]). Each graph corresponds to a different minicircle size out 
of L = 1.2, 2.9, 5.8, and 11.9 x (persistence length). 
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a = (1.00 



a = 0.05 



a = 0.13 



Ljlp « 1.2 






FIG. 7: Snapshots of model DNA minicircles at three different lengths and three different excess 
linking numbers densities (cj). Note that some of the bonds are omitted in the figure for visual 
clarity. 

AL/c among writhe and twist depends on the stiffnesses associated with the two quantities 
(i.e., bending and twisting rigidities) which are complex functions of the intra- and inter- 
strand interactions between the superatoms. Let 

ALA; Lk . 



(9) 



where positive and negative values of o correspond to overtwisted and undertwisted circular 
DNA chains, respectively. Note that, chains of different size are under similar local torsional 
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stress if their a values are identical. Nevertheless, we show below that the response of the 
DNA to cr 7^ depends strongly on the length. 

As discussed above, the circular model DNA chain carries a finite writhe {Wtq) and twist 
{Twq) in the relaxed state [Lk = Lko) due to the corkscrew motion of the centerline. Since 
we are interested in measuring the change in twist and writhe as a function of cr, we define 
the average deviation in twist as 

(ATw) 
Lko 

and the average deviation in writhe as 

(Am-) 

Lko 

where the averages are taken over equilibrium snapshots of the system. 

In Fig. [6], (AWr) / Lko and (ATw) / Lko are given as a function of a. Four different chain 
lengths are considered, L/lp ^ 1.2, 2.9, 5.8, and 11.9 and the results are shown in separate 
graphs from top- left to bottom-right, respectively. 

Since Lk is an integer, the values of a realizable for a fixed chain length L are discrete. 
In order to overcome this constraint, we combined data obtained from chains with lengths 
varying upto ±5% of the chosen L/lp while sampling cr g [—0.20, 0.20] by changing Lk. For 
example, L/lp ^ 1.2 regime was sampled with {Lk = 7,8, ... , 12} (^{L = 105, 106, . . . , 117}. 
Above variability in length does not give rise to a significant error in the average writhe and 
twist of long chains. For short chains the effect is more pronounced and fluctuations are 
observed in the data. Nevertheless, a general trend that varies with the chain length is clear 
in Fig. [6] and will be discussed next. 

The average writhe and twist of the circular chains display a nonmonotonic dependence on 
a for all four cases considered. Let us first focus on the shortest chain regime {L/lp ^ 1.2) 
shown in the top- left graph of Fig. [6] and the overtwisting scenario with a > 0. There 
exist two qualitatively distinct modes of torsional response which are separated by a sharp 
transition at cr+ ^ 0.07. For small deviations from Lko with a < a^, the extra linking number 
is completely absorbed by the change in twist. In this window the minicircle essentially 
remains planar (Fig. [71 top-left), with a slightly negative slope in writhe which is probably 
associated with the aforementioned nonzero writhe density of the relaxed chain. 

At the transition, the twist drops sharply and the writhe increases accordingly, absorb- 
ing both the additional linking number and the reduction in twist. The jump in writhe 
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(Tw) - Twq 
Lko 



{Wr) - Wro 
Lko 



(10) 



(11) 



is manifested as an out-of-plane deformation of the DNA minicircle (Fig. [71 top-middle). 
This phenomenon corresponds to the supercoihng transition of the minicircle which can be 
understood through the buckling instability of the planar state in a circular elastic ribbon 
model. In fact, the present data can be used to extract the twist persistence length (= K), 
of our model DNA through the relation 32 1 



K-Lko 

for an elastic ribbon with a symmetric torsional response. Substituting a, 



(12) 



0.07 yields 

K/lp ^ 2.5 , consistent with earlier analysis [33| of the single-molecule experiments. {4] 
Further increase in writhe leads to a complete figure-eight shape in the DNA chain as seen 
in Fig. [7] top-right. Beyond the buckling point, the linking number is almost equally shared 
by writhe and twist. 

A qualitatively similar behavior is observed upon undertwisting in the interval > a > 
—0.1. The buckling transition upon undertwisting takes place at a" ~ —0.09. The asymme- 
try I (7+ 1 < \a^\ suggests that undertwisting is easier than ovetwisting the DNA chain. Such 
nonlinear response in torsional stiffness has already been reported in experiments 34, l35|] and 
full-atomistic cornputer simulations 36||. A recent analytical model for DNA minicircles by 
Liverpool et al. 16| also predicts that DNA minicircles favor supercoihng to denaturation 
in the weakly nonlinear regime, consistent with Fig. [61 A possible origin of the nonlinearity 
is the presence of twist-bend coupling in chiral molecules where the symmetry under 180" 
rotation around the helix axis is broken (i.e., molecules with major and minor grooves) as 
argued by Marko&Siggia. [s^] Note that, we do not see local base-pair hydrogen-bond break- 
ing events ISj in this regime. Since our training data set does not include any information 
regarding the melting of the double helical structure, it is not surprising that the coarse 
grained model also does not display melting. Incorporation of full-atomistic single-strand 
DNA simulation data into the model could probably suffice to capture such behavior, which 
is planned as future work. 

When we compare the L/lp ^ 1.2 case to longer chains, we observe two significant 
differences. First, the buckling transition takes place at smaller lad values, consistent with 
Guitter and Leibler 32[, and practically disappearing for L/lp > 6.0. At constant torsional 
stress density (e.g., a = 0.05), minicircles with increasing length accommodate a higher 
fraction of ALk in writhe, with typical configurations shown in the middle column of Fig. [3 
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Second, for L/lp 2.9, 5.8, and 11.9 beyond the buckling point, the majority of the 
excess hnking number is absorbed by the writhe, unhke the equal distribution we have seen 
for L/lp 1.2. For these longer chains, {AWr) / {ATw) decreases with increasing |cr|, which 
is another manifestation of the nonlinear twist rigidity. For fixed cr, (AW^r) / (ATw) increases 
monotonically with length, possibly approaching a finite asymptotic value. 



D. Summary and Conclusions 

We presented a coarse-grained model which is designed for studying the equilibrium 
structural properties of 10^-10^ bps long DNA minicircles. Proper thermodynamic averaging 
of global structural features requires microsecond simulations. Therefore we chose a minimal, 
two-bead representation of the sugar-phosphate and the base in the basic repeat unit. This 
approach may be used to address single-strand chirality and base-flipping which are not 
accessible to single-bead models. The model parameters were extracted from full-atomistic 
molecular dynamics simulations of DNA oligomers via Boltzmann inversion. Even at this 
level of simplicity, used coarse-graining methodology yields a faithful representation of the 
directionality, helicity, major and minor grooves and similar local characteristic features of 
DNA. For the sake of simplicity and efficiency, base-pair specificity and explicit electrostatic 
interaction have been ignored here, althought it is straightforward to incorporate these 
into the model within the present approach. Denaturation upon underwinding should be 
observable after adding base specificity and including full-atomistic simulations of single- 
strand DNA in the Boltzmann inversion step, which will be considered in a future extension. 

Using our model, we performed a systematic molecular dynamics study of supercoil forma- 
tion in DNA minicircles. In particular, we measured the twist/writhe partitioning expressed 
in Eq. ([H]) as a function of the chain length (L) and excess linking number density (a). 
We observed a supercoiling (buckling) transition associated with the off-plane deformation 



instability of the minicirc 
a simple elastic model [16 



es 



br > \af{L)\, as predicted by analytical calculations on 
32| and recent full-atomistic simulations, [l^ The transition is 



marked by a sudden increase in the writhe, while {ad decreases with L and practically disap- 
pears beyond L/lp ^ 6 (cr^ ~ L^^ for elastic models). In the planar regime with \a\ < |cr^|, 
the excess linking number is essentially stored in the twist. 

Our results suggest that, beyond the supercoiling transition, the fraction of the linking 
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number absorbed as twist and writhe is also nontrivially dependent on chain length. Chains 
of the order of a persistence length carry approximately equal amounts of twist and writhe, 
while longer chains accommodate an increasing fraction of the excess linking number as 
writhe. The dependence of (ATw) / (AWr) ratio on a remains nonlinear also for larger 
chains. At fixed a, we observe that the above ratio increases with L and possibly reaches 
a finite asymptotic value in the limit L — oo, as this limit is accurately represented by 
harmonic elastic energy terms for both twist and writhe. On the other hand, the behavior 
at fixed a/L, which may apply to the bound portions of a circular DNA chain near the 
melting temperature, is unclear and appears to be an interesting problem. 
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